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Abstract 



I J* In this work we derive and analyze variational integrators of higher order for the structure-preserving 

^H simulation of mechanical systems. The construction is based on a space of polynomials together with 

(-H Gauss and Lobatto quadrature rules to approximate the relevant integrals in the variational principle. 

■4-J The use of higher order schemes increases the accuracy of the discrete solution and thereby decrease 

2 the computational cost while the preservation properties of the scheme are still guaranteed. The order 

f^ of accuracy of the resulting variational integrators are determined analytically and it is shown which 

' ' combination of space of polynomials and quadrature rules provide optimal convergence rates. For partic- 

, ular integrators the order can be increased compared to the Galerkin variational integrators previously 

K^ introduced in [MWOl]. Furthermore, linear stability properties, time reversibility, structure-preserving 

QQ properties as well as efficiency for the constructed variational integrators are investigated and demon- 

^v strated by numerical examples, discrete variational mechanics; optimal order rates; symplectic methods; 

^f^ variational integrators 



^ 1 Introduction 

^ — I 

. 5^ During the last years the development of geometric numerical integrators has been of high interest in numeri- 

\^ cal integration theory. Geometric integrators are structure-peserving integrators with the goal to capture the 

5—1 dynamical system's behavior in a most realistic way ([MWOl, HLW02, Rci94]). Using structure-preserving 

methods for the simulation of mechanical systems, specific properties of the underlying system are handed 
down to the numerical solution, for example, the energy of a conservative system shows no numerical drift or 
the first integrals induced by symmetries are preserved exactly. One particular class of structure-preserving 
integrators is the class of variational integrators, introduced in [MWOl, Sur90], and which has been further 
developed and extended to different systems and applications during the last years. Variational integrators 
([MWOl]) are based on a discrete variational formulation of the underlying system, e.g. based on a discrete 
version of Hamilton's principle for conservative mechanical systems. The resulting integrators are symplectic 
and momentum-preserving and have an excellent long-time energy behavior. 

By choosing different variational formulations (e.g. Hamilton, Lagrange-d'Alembert, Hamilton-Pontryagin, 
etc.), variational integrators have been developed for a large class of problems: These involve classical 
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conservative mechanical systems (for an overview see [LMOW04a, LMOW04b]), forced and controlled sys- 
tems ([OBJMll, KMOWOO]) , constrained systems (holonomic ([LMO08, LOBMOIO]) and nonholonomic 
systems ([KMSIO, CMOl])), nonsmooth systems ([FMOW03]), stochastic systems ([BRO08]), multiscalc sys- 
tems ([TOMIO, L0B13, SG09]) and Lagrangian PDE systems ([LMOW03, MPS98]). The appHcability of 
variational integrators is not restricted to mechanical systems. In [OBTC+13] variational integrators have 
been developed for the structure-preserving simulation of electric circuits. The generation of different meth- 
ods, e.g. Lie group variational integrators ([LLM07]), as well as extensions to different system classes are 
recent research topics in variational integration theory. 

Of special interest is the construction of higher order symplectic integrators: To ensure moderate computa- 
tional costs for long-time simulations, typically first or second order integrators are used. However, many 
applications, in particular problems in space mission design, demand more accurate discretization schemes. 
There are mainly three different ways of constructing symplectic integrators of higher order (cf. [LR04]): 
(i) By applying composition methods higher order symplectic schemes can be constructed in a system- 
atic way based on a splitting of the Hamiltonian into explicitly solvable subproblems (for an overview see 
[MQ02, Yos90]). (ii) For (partitioned) Runge-Kutta methods there is a well-developed order theory which 
can be used to identify higher order symplectic Runge-Kutta methods (order conditions on the coefRcicnts 
have first been introduced by [SS88, Las88, Sur89, Sun93]). However, the identification of the coefficients 
for a symplectic Runge-Kutta scheme of desired order is not trivial and quite involved, (iii) In contrast, 
generating functions (see e.g. [Arn99]) can be constructed that automatically guarantees the symplecticity 
of the associated numerical method. Since variational integrators rely on the approximation of the action, a 
generating function of first kind (cf. [MWOl]), we focus on the latter approach for the construction of higher 
order methods. 

Of particular interest are Galerkin variational integrators which have been already studied in e.g. [MWOl, 
Leo04, LSlla, HL12]. They rely on the approximation of the action based on a choice of a finite-dimensional 
function space and a numerical quadrature formula. In [LSlla] Galerkin and shooting-based constructions 
for discrete Lagrangian are presented. Rather than choosing a infinite-dimensional function space, the 
shooting-based construction depends on a choice of a numerical quadrature formula together with a one-step 
method. For example, the prolongation-collocation variational integrators based on Hermite polynomials 
and introduced in [LSllb] fall into this category. In [HL12] a convergence analysis for Galerkin type varia- 
tional integrators is established showing that under suitable assumptions the integrators inherit the order of 
convergence given by the finite-dimensional approximation space used in the construction. More detailed, it 
is shown that a Galerkin variational integrator based on an s -I- 1-dimensional function space (e.g. the space 
of polynomials of degree s) and a quadrature rule of order s has convergence oder s. 

In this contribution we show that the convergence order of the variational integrator can even be increased 
if higher order quadrature rules are used. We focus on a particular case of Galerkin variational integrators: 
As finite-dimensional function space we choose the space of polynomials of degree s. Furthermore, as 
quadrature rules we focus on the Gauss and Lobatto quadrature formula. However, in contrast to [MWOl] 
we do not restrict the number r of quadrature points being equal to the polynomial order s. By using 
techniques from variational error analysis we investigate which combination of space of polynomials and 
quadrature rules provide optimal convergence rates. In particular, we prove that the order of the higher order 
variational integrator constructed by a polynomial of degree s and a quadrature rule of order u is min (2s, u) 
(Theorem 4.6). Thus, the integrator order can be increased to 2s for sufficient accurate quadrature rules. 
Furthermore, we investigate preservation properties, time reversibility and linear stability of the constructed 
integrators and perform a numerical analysis regarding efficiency versus accuracy by means of numerical 
examples (see also [Saal2]). 



Outline In Section 2 we recall the basic definitions and concepts of variational mechanics and variational 
integrators. In Section 3 first the higher order integrators are constructed following the Galerkin approach 
introduced in [MWOl] and preservation properties, such as symplecticity and preservation of momentum 
maps are discussed (Section 3.4). It is shown under which conditions the constructed integrators are time- 



reversible (Section 3.5). Furthermore, a linear stability analysis is performed for specific examples showing 
in which region the constructed higher order variational integrators are asymptotically stable (in the sense 
that the growth of the solution is asymptotically bounded (cf. [LR04]), Section 3.6). The main contribution 
of this work is given in Section 4. The concept of variational error analysis is introduced and applied to prove 
the optimal order of accuracy of the higher order variational integrators. In Section 5 the analytic results 
are demonstrated by two numerical examples, the harmonic oscillator and the Kepler problem. A numerical 
analysis regarding efficiency versus accuracy is performed for two different classes of variational integrators. 
Finally, we conclude with a summary of the results and an outlook for future work in Section 6. 



2 Variational mechanics 

2.1 Hamilton's principle and Euler-Lagrange equations 

Consider a mechanical system defined on the n-dimensional configuration manifold Q with corresponding 
tangent bundle TQ and cotangent bundle T*Q. Let q[t) G Q and q{t) € T^^^f-^Q, t € [0,T] be the time- 
dependent configuration and velocity of the system. 

The action & : C^([0,T],Q) — > M of a mechanical system is defined as the time integral of the Lagrangian, 
i.e., 

&{q)= f L{q{t)A{t))dt (2.1) 

where the C^-Lagrangian L : TQ — > M consists of kinetic minus potential energy. Hamilton's principle seeks 
curves q e C^([0,T],(5) with fixed initial value q(0) and fixed final value q{T) satisfying 



6&{q) = (2.2) 



for all variations 5q G TqC'^([Q,T],Q). This leads to the Euler-Lagrange equati 



ions 



which are second-order differential equations describing the dynamics for conservative systems. 



2.2 Discrete Hamilton's principle and discrete Euler-Lagrange equations 

The concept of variational integrators is based on a discretization of the variational principle (2.2). Consider 
a time grid At = {tk ~ kh\k — Q, . . . ^ iV}, Nh — T, where iV is a positive integer and h the step size. We 
replace the configuration q{t) by a discrete curve qd — {qk}k=() with q^ = qd{tk) as approximations to q{tk)- 
We define a discrete Lagrangian L^ : Q x Q — >■ M 

Ldiqk,qk+l) ~ / '^' L{q{t),q{t)) dt (2.4) 

Jtk 

that approximates the action on [i^, t^+i] based on two neighboring discrete configurations q^ and qu+i- The 
discrete action <3d : Q^^^ ^ M is defined as 

N-l 

&d{qd) = X! ^d{qk,qk+i)- 

k=0 

The discrete Hamilton principle is formulated by finding stationary points of the discrete action given by 

dediqd) = (2.5) 



with Sqo = Sqn = 0. This gives the discrete Euler- Lagrange equations (DEL) 

DiLd(qk,qk+i) + -D2^d(gfe-i, 9fe) = (2.6) 

for k — 1, . . . ,iV — 1 and with Di being the derivative w.r.t. the i-th argument. Equation (2.6) provides 
a discrete iteration scheme for (2.3) that determines qk+i for given qk-i and qk- It is also known as the 
discrete Lagrangian map Fl^ : QxQ -)■ QxQ, given by FL^iqk-i,qk) = iqk,qk+i) and {qk-i,qk), {qk,qk+i) 
satisfy (2.6). The discrete iteration schemes derived by a discrete variational principle are called variational 
integrators and are well-known to be symplectic and momentum-preserving and exhibit excellent long-time 
energy behavior (cf. Section 3.4). 

The discrete Legendre transforms ¥ Ld : Q x Q ^?' T*Q provide discrete expressions for the conjugate 
momenta by 

¥~Ld ■■ {qk,qk+i) -^ i(lk,Pk) = i(lk,-DiL{qk,qk+i)) and 
V'^Ld : (gfe-i,(?fe) -> {qk,pt) = {qk:D2L{qk-i,qk))- 

Note that (2.6) can be equivalently written as p^ — p^ and the discrete Hamiltonian map Fj^^ : T*Q ^ T*Q 
defined by 

Fl^ ■■ {qk,Pk) -^ {qk+i,Pk+i) = F±id o Fl^ o (F±Ld)-i((jfc,pfe). 

is equivalent to the discrete Lagrangian map. 

In [MWOl] a relation between the error of the discrete Hamiltonian map compared to the analytic Hamil- 
tonian flow and the error of the discrete Lagrangian compared to the analytic action (see (2.4)) is given, 
whereas the latter approach is denoted as variational error analysis. In Section 4 we use this result to prove 
the orders of accuracy of the variational integrators constructed in this work. 



3 Higher order Variational Integrators 

3.1 Quadrature rules 

For the approximation of the action arbitrary quadrature rules can be used. Let us briefly repeat the basic 
definitions (cf. e.g. [Hil87, SK06]). In the following, the quadrature formula /(/) = J2l=i'^ifi^i) given by 
the weights Wi and the quadrature points Xi, i = 1, . . . , r, for the approximation of /(/) = / f{x) dx is 
denoted by quadrature formula {wi,Xi)l^i w.r.t. the interval [a,b]. 

Definition 3.1 (Accuracy of quadrature formula). The quadrature formula I{f) = ^^^iWifixi) given 
by the weights Wi and the quadrature points Xi for the approximation of L(f) = J f{x)dx has a degree of 
accuracy of k Cz Nq if it integrates all polynomials of degree smaller or equal to k exactly and k is the maximal 
number satisfying this property. 

In particular, we have, 

Hp) ~ /^ Wip{xi) — I p{x) dx = L{p) for all p G H', / < /c, 
i=i "'" 

with n' being the space of polynomials of degree I. In this work, we use Gauss and Lobatto quadrature rules. 
The Gauss quadrature has maximal degree of acuracy 2r — 1 for a fixed number r of quadrature points, while 
the Lobatto quadrature has degree of accuracy 2r — 3 for the same number of quadrature points. Quadrature 
points and weights of the Gauss- and the Lobatto quadrature formula w.r.t. the interval [—1,1] are listed in 
Table 1 (cf. [Zwi03]). The order of a quadrature formula is defined as follows. 



Table 1: Weights and quadrature points of the Gauss- and the Lobatto quadrature formula w.r.t. the interval 
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Gauss quadrature Lobatto quadrature 
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Definition 3.2 (Order of quadrature formula). The quadrature formula {wi,Xi)l^i w.r.t. the interval [0,1] 



«-i 



is of order u if for all f G Cli ^, we have 

h 



f 

/ /(x)dx = /iVwJ(a;,/i)+C'(/i"+i) for h -^ 0. 
•^0 ,=1 



The relation between order and accuracy of a quadrature formula is given by the following theorem. 
Tiieorem 3.3. A quadrature formula is of order u if its degree of accuracy is u — 1. 



Proof. Let / € ^Tqu] ^"^^ ^^* ^^^ quadrature formula {wi,Xi)\^i w.r.t. the interval [0,1] be of degree of 
accuracy u — 1. Note that the quadrature formula can be transformed to an arbitrary interval [a,b] with 
weights Wi = j^ and quadrature points Xi — ^^E^ ■ The Taylor expansion of / in x = is given by 

^2 ^M— 1 

/(x) =. /(0) + x/'(0) + -/"(0) + ... + ^^^-^/(«-i)(0) + O(x") 

:= T"-\/(x)+0(x"), 

where we denote by T"^^/(x) the Taylor polynomial of / in up to order u — 1. T"^^/(x) is a polynomial 
of degree u — 1 and thus can be integrated exactly by the applied quadrature formula. It follows 

^f{x)dx = ^'(^/(0)+x/'(0) + |[/"(0) + ... + ^^^/("-i)(0)+O(x"))dx 

h j-l 



T"- V(a;) dx + C'(/i"+^) = h / T'^-'f{xh) dx + 0{h'^+') 
/o Jo 

r r 

^ hY, w,T^-^f{x,h) + 0{h"+^) ^hJ2 mifix^h) + 0{{x,hY)) + 0(/i"+i) 

r 



3.2 Approximation of the action integral 

The approximation of the action integral consists of two approximation steps: the approximation of the space 
of trajectories and the approximation of the integral of the Lagrangian by appropriate quadrature rules. 

We approximate the space of trajectories C{[0,h],Q) = {q : [0,h] ^ Q \ g(0) — qa,q{h) — qb} by a finite- 
dimensional approximation C*([0, h],Q) C C([0, h],Q) of the trajectory space given by 

C'{[0,h],Q) = {qeC{[0,h],Q)\qen^}, 

with n'* being the space of polynomials of degree s. Given s + 1 control points = do < di < ■ ■ ■ < dg-i < 
ds = 1 and s + 1 configurations go = {Qq, QotQq, ■ ■ ■ , QcT ^lo) ^ Q''^^ with (7° = qa and ^q = qb, the degree s 
polynomial qd{t; qo, h) which passes through each q^ at time d^h, that is, qd{d^h) ^ q^ ior v — 0, . . . , s, is 
uniquely defined. 

With the Lagrange polynomial Z^.s : [0, 1] — > M 

we obtain qd(t; qg, h) with t G [0, h] as 

qd{t;qo,h) ^^q^l^^ \h)' 
With 

for all z = 0, . . . , s, we have qd{dih] qQ,h) = ^q. The derivative of qd{t', qo, h) w.r.t. t provides an approximation 

of q on [0, h] as 



1 -ft 

i{t;qa,h) = -r^qoL.s I t 



To approximate the trajectory q : [0,T] — )■ Q, we divide the time interval [0,r] in iV = T/h sub intervals of 
length h as 

AT-l 

[o,r]= IJ [kh,{k + i)h]. 

fe=0 

On all sub intervals we approximate q : [0,T] — > Q piecewise by the polynomials qd,k ■ [kh, {k + l)h] — > Q 
defined by 

grf,fc(i; qk,h) ■= qd{t - kh; qk,h) t € [kh, (fc + l)h] 

with qk = {q1, q]., ■■■, ql~ , qf.) and fc = 0, . . . , A^ — 1. To obtain a continuous approximation on [0, T], we set 
ql = qt+,ioTa\lk = 0,...,N-2. 

For the approximation of the action integral (2.1) we replace the curves q{t) and q{t) by the piecewise 
polynomials qd,k{t', qk, h) and qd,k{t', qk,h), k = 0, . . . ,N, and approximate 

k+l)h 

L{qd,k{t; qk, h), qd.k{t; qk, h)) dt (3.1) 

kh 

on each time interval [kh, {k + l)h], k — Q, . . . ,N — 1, by choosing a numerical quadrature rule {hi, Ci)\^i 
w.r.t. the time interval [0,1] with quadrature points ci £ [0,1] and weights bi, i = 1, . . . ,r. The choice 
of quadrature rule should be adapted to the desired order of accuracy of the integrator since the or- 
der of the quadrature rule provides an upper bound for the order accuracy of the variational integrator 



(cf. Theroem 4.3). By applying the quadrature rule (bi,Ci)l^i to the integral (3.1) we define the discrete 
Lagrangian Ld,k as 

r 

Ld,k = Ldiqk = iql,—,qt:),h) = h^biL{qd,k{cih + kh; qk, h), qd,k{cih + kh; qk, h)) 

i=\ 
r 

= h'^biL{qd{cih; qk,h),qd{cih; qk, h)) 

i=l 

which provides an approximation of the action on the interval [kh, {k + l)/i] as Ld^k ~ /i.i L{q{t) , q{t)) dt. 

Note that the discrete Lagrangian depends on s + 1 configurations qk and the step size h. In the following, 
we write Ld{qk) for Ld{qk, h). Finally, we define the discrete action sum over the entire trajectory to be 

JV-l 

&d{qo,---,qN-i)^^Ld{qk) (3.2) 

fc=0 

which is an approximation of the action sum on [0, T] as &d{qo, ■ ■ ■ , qN-i) ~ ®{q)- 

3.3 Discrete Hamilton's principle 

For the discrete action &d defined in (3.2) we can apply discrete Hamilton's principle as described in Sec- 
tion 2.2. Since we want to determine discrete approximations of curves for which the discrete action is 
stationary, the derivatives of the action w.r.t. q^ have to vanish for all fc = 0, . . . , A'^ — 1 and z^ = 0, . . . , s. 
This leads for fc = 0, . . . , iV — 1 and i' — 1, ... ,s — 1 to 

= _(,„,...,g^_,)^_(g,) 



h'^bi ( -^{cih;qk)L,s{ci) + -wr{cih\qk)-iu,s{ci) j 



Note that we use the short notation i-{cih;qk) = ^{qd{cih;qk,h),qd(cih;qk,h)) and that we have jprr = 
dqdc^h^gk.h) ^ /^^(c^). The analog holds for the other two terms. With ql_-^ = q° for all /c == 1, . . . , TV - 1 
we obtain for z^ = and ly = s 



d&d , . d&d, . dLd{qk-i) , dLd{qk) 

= g^^{qo,-,qN-l)-^{qo,...,qN-l)-^^^-- + ^^ 

= h'^bi i —{cth;qk-i)ls,s{ci) + -7y{cih;qk-i)-is,s{ci) 

dL 
dq 



'h'^bi I -T-{cih;qk)lo,s{ci) + -Trricih;qk)-io,s{ci 



With the notation DiLd{q1, . . . ,q^) ■— ''^ _ i' we obtain the discrete Euler- Lagrange equations 

"Ik 

D,+iLd{qt^„ ..., ql_,) + D,Ld{ql, ...,ql) = 0, (3.3) 



Aidfe°,...,9D = Vz = 2,...,s, (3.4) 

for fc = 1, . . . , A^ — 1. A sequence {qk}k=o = iilk^ ■■■i1k)}k=o ^^^^ satisfies (3.3)-(3.4) and the transition 
condition is a solution of the discrete Euler- Lagrange equations. We denote the left hand side of (3.3) and 
(3.4) by DDEhLdiqk-i,qk) such that we have 

{DDELLd{qk-i,qk))i == Ds+iLd{qk-i:--;q'l^i) + DiLd{q°,...,ql) 
{DuEhLd{qk-i,qk))2 = ^2id(g°, •■•, gfc) 

{DuELLd{qk-i:qk))s = DsLd{ql,...,ql). 

As in [MWOl] we can introduce the standard discrete Lagrangian that depends only on two configurations 
as 

Ldiql^ql) = Ldiqk), 

where q^, . . . , g^~ are implicitly determined by satisfying the internal stage equations (3.4). Alternatively, 
one can characterize the discrete Lagrangian in the following way (see [MWOl]), 

r 

Ld{ql,q'l)= ext h^biL{qd{cih;qk)),qd{cih;qk)), 

!^e{i,...,s-i} 

meaning that (s— 1) configurations q^, . . . , q^~ are determined by extremizing the discrete Lagrangian. The 
Lagrangian Ld{q^,ql^) provides the same iteration scheme as the discrete Lagrangian Ld{qk)- 

Let (7fe — (g^,...,5|) and let a{qk,qk+i) = <7fe+i be the translation operator and TT{qk,qk+i) = qk the 
projection operator. The discrete Lagrangian evolution operator X^^ : Q^^^ — >■ Q^^^ x Q'^'^^ satisfies 

XLa{qk-i) == {qk-i,qk) with g^_i = ql, 

T^ ° XLaiqk-i) ^ qk-i and -DoEL-^dO A:L^((?fe-i) = 0. 
The discrete Lagrangian map F^^ : Q^^^ — > Q*^^ is defined by 

FhAqk-i) =aoXLa{qk-i) = qk 

and generates the sequence of configurations that is denoted as the solution of the Euler-Lagrange equations. 
The discrete Legendre transforms ¥^Ld : Q x Q —>■ T* Q are defined as 

¥-Ld{qlqi) = {qlpln = {ql-D,Ld{qlqi)), 

F+irf(gLi,g|_i) = iqlpl+) = iqlD,+,Ld{qt_„ql_,)). 

From the discrete Euler-Lagrange equations it follows that along the solution of the discrete Euler-Lagrange 
equations we have 

Pk •- Pk - Pk ■ 

Note that the discrete Lagrangian flow is well-defined, if the discrete Lagrangian is regular, i.e., if the discrete 
Legendre transforms are local isomorphisms what is assumed in the following (see also Assumption 4.4). As 
shown in Fig. 3.1, the discrete Hamiltonian map F^^ : T* Q -^ T* Q is given by 

FL,=¥^LdoFL,o{¥^Ld)-\ 



With the proposed method different variational integrators can be constructed. We use the following no- 
tation: (PsNrQu) is an integrator constructed as described above with a polyomial of degree s with sym- 
metrically distributed control points (di),f^o ^^'^ ^ quadrature formula of order u with r quadrature points. 
Note that u depends on r. If explicitly given, we denote by three letters the quadrature rule in use, i.e. Lob 
for Lobatto quadrature {u = 2r — 2) and Gau for Gauss quadrature {u = 2r). 



(9§,--->9g)h 



Fl, 



-^iql...,qf) 




(qIp'o) h 



Fr 



-^ [qWi) h 



Ft, 



-^ {qlpD 



Figure 3.1: Correspondence between the discrete Lagrangian and the discrete Hamiltonian map. 



Example 3.4. (i) The integrator {PlNlQ2Gau) is based on a polynomial with control points do = 0, 
di = 1 and the quadrature approximation 

L{q{t),q{t))dt^Ld{{qlq°),h) ^ hL{qd{h/2) Ad{h/2)) 



= HL 



%" + '?? 9?-%° 



and thus results in the midpoint rule discrete Lagrangian. 

(ii) If the trapezoidal rule is applied as quadrature formula, we obtain the variational integrator {PlN2Q2Lob) 
with discrete Lagrangian 

La{{qlql),h) = ^L{qM,qdm + \L{qa{h)Ad{h)) 






h . ( q\ - ql 



r^ Ui 



2 V h 

which is the discrete Lagrangian of the Stormer- Verlet method with factor ^ . 

(Hi) For a second order polynomial (s — 2) with control points do = 0, di = ^ and d2 = 1 and by 
applying Simpson's rule, which is a Lobatto quadrature of order four, we obtain the variational integrator 
(P2N3Q4:Lob) with discrete Lagrangian given by 

Ldiiql 9o, 9o), h) = ^i(gd(0), 9rf(0)) + ^L{qd {h/2) , q, {h/2)) + \L{qa{h) , qd{h)) 

DO 

h^ f ^ -3(?0 + ^ql -q^\ 2hr ^ ql-q°„\ hr ^ q^o- H + ^Qo 



- 6^(^o 



3-^'^- h 



+ qL ( ql 



Remark 3.5 (Implementation). For given configurations q^,. . . ,q^ ,qi we compute the unknown configu- 
ration 52 by performing one step of the discrete Lagrangian evolution 

{q'„...,q^,-\q',)^{ql...,qr\q',). 

To this end, the system of s nonlinear equations (3. 3)- (3. 4) is solved for the s unknowns {ql, ..., g*~ , (?2)- For 
the numerical solution, a Newton method can be used. For given initial configuration q{0) and momentum 
p{0), in the first .step, (3.3) is replaced by the discrete Legendre transform 

p(0) = -i?iLd((jO,...,9^) 

and the system of equations is solved for {q^, . . . , q\). For the quadrature rules considered here, (3.3)-(3.4) 
typically provide an implicit scheme for nonlinear systems that has to be solved by an iterative solver all at 
once. 



Remark 3.6 (Galerkin methods). For the Galerkin variational integrators as introduced in [MWOl], Sec- 
tion 2.2.6, the number of quadrature points of the quadrature formula (bi,Ci)l^i is fixed to r = s, where s 
is the degree of the polynomial qd. In our notation that means that only the methods (PsNsQu) are inves- 
tigated, which are shown to be equivalent to partitioned sym,plectic Runge-Kutta methods. In particular, it 
is pointed out that the integrator {PsN sQ2sGau) , which uses the Gauss quadrature formula, corresponds to 
the collocation Gauss-Legendre rule, whereas {PsNsQ2s — 2Lob) yields the standard Lobatto IIIa-IIIB par- 
titioned Runge-Kutta method. For both methods the order is determined by the order of quadrature rule, i.e., 
{PsNsQ2sGau) is of order 2s and {PsNsQ2s - 2Lob) is of order 2s -2 (cf. [HLW02]). Although the Gauss 
quadrature formula leads to higher order schemes, in particular for stiff systems the choice of a quadrature 
rule involving Cg = \ leads to better numerical performance (cf. [HLW02, MWOl]). If, in addition, one 
wishes to use a symmetric quadrature rule, i.e., ci — 0, the Leg endre- Lobatto quadrature rule provides the 
highest possible order. A detailed order analysis is performed in Section 4- 



3.4 Preservation properties of variational Integrators 

As already mentioned in Section 2 variational integrators are structure-preserving, in particular they are 
symplectic and momentum-preserving. In this section we briefly repeat the notion of symplecticity and first 
integrals and their discrete counterparts. 



3.4.1 Symplecticity and energy behavior 



In the case of conservative systems (as considered here), the flow on T*Q of the Euler-Lagrange equations 
preserves the canonical symplectic form n = dq^ A dp^ = A9 of the Hamiltonian system, where 9 = p^dq^ is 
the canonical one-form. It is well known (cf. e.g. [MWOl]) that variational integrators are symplectic, that 
is the same property holds for the discrete flow of the discrete Euler-Lagrange equations. As a consequence, 
the canonical discrete symplectic form f2 = dqQ A dp^ is exactly preserved for the discrete solution which is 
in particular true for the variational integrators presented in this work. 

By using techniques from backward error analysis it is shown (cf. e.g. [HLW02]) that symplectic integrators 
have excellent energy properties, meaning that for long-time integrations there is no artificial energy growth 
or decay due to numerical errors. This property is demonstrated numerically in Section 5 and illustrates 
the advantage of variational integrators over e.g. nonsymplectic Runge-Kutta integrators in particular for 
long-time simulations. 



3.4.2 Preservation of first integrals 

The Noether theorem provides first integrals of the Euler-Lagrange equations which are also called momentum 
maps. The discrete Noether theorem states that these invariants are also preserved for the discrete solution. 
The following two theorems are taken from [HLW02] . 

Theorem 3.7 (Noether theorem). Let L{q,q) be a regular Lagrangian. Suppose G — {g^ : v e M} is a 
one-parameter group of transformations (gy o g^, = g^+w) which leaves the Lagrangian invariant such that 

L{9v{q),9'M)q) = L{qA) V« G M V(q, (?) G TQ. 

Let a[q) = -^gv{q)\v=o be defined as the vector field with flow gv{g)- Then 

I{q,p)=p^a[q) (3.5) 

is a first integral of the Euler-Lagrange equations. 
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The discrete analog of the Noether theorem is stated as follows. 

Theorem 3.8 (Discrete Noether theorem). Suppose the one-parameter group of transformations G ~ {g^ : 
V G M} leaves the discrete Lagrangian Ld{qQ,qi) invariant, that means 

Then (3.5) is an invariant of the discrete Hamiltonian map F^^, i.e., 

IoFL,{qlpl)^I{qlpl). 

Proofs of Theorems 3.8 and 3.8 can be found in [HLW02]. 

Note that the invariant of the discrete Hamiltonian map only equals the first integral of the Euler-Lagrange 
equations if the discrete Lagrangian Ld inherits the invariance of the Lagrangian L. In the following we show 
under which condition this invariance is inherited. 

Definition 3.9 (Equivariance). Let {g^ : v G M} be a one-parameter group of transformation. The interpo- 
lation polynomial qd of {PsNsQu) is equivariant w.r.t. {gy : v G M} if 

9v{qd{t;{qo,...,qo))) = qd{t;{gv{qo),---,9v{qo))), (3.6) 

g'y{qd{t; {q"o, . . . , q'o)))q{t; {ql . . . ,gg)) = gd(i; (5.(90°), • ■ ■,9v{q'o)))- (3.7) 

Note that (3.7) follows from (3.6) by applying the chain rule. 

Theorem 3.10 (Invariance of discrete Lagrangian). Suppose that the interpolation polynomial qd of {PsNrQu) 
is equivariant and that the regular Lagrangian L is invariant w.r.t. the one-parameter group G — {g^, : v G M}. 
Then the discrete Lagrangian Ld of (PsNrQu) is invariant w.r.t. G. 

Proof. Let go = (^o: • • • :9o) with q^ = q^ and g„ • go = (g«(9o), ■■■,9v{qr^),9v{qi))- Let {hi,Ci)\^^ be the 
quadrature formula that corresponds to (PsNrQu). With the invariance of L and the equivariance of qd we 
have that 

r 

Ld(9v{qo),9v(qi)) = ext h^biL(qd(cih;gy ■ qa),qd(cih;gy ■ qo)) 

iye{l,...,s-l} 

r 

= ext h^biL(g^{qd{cih;qo)),gy(qd(cih;qQ))qd(cih;qQ)) 

iye{l,...,s-l} 

r 

= ext h^biL(qd(cih;qo)),qd(cih;qo)) 

u£{l,....,s-l} 

- Ld(q"o,q"i)- 



A general form of the statement of Theorem 3.10 for Galerkin Lie group variational integrators can also be 
found in [LSlla]. 

Remark 3.11. From Theorem 3.10 it follows that for linear group transformations the discrete Lagrangian 
Ld of (PsNrQu) inherits the invariance of the Lagrangian L (cf. [MWOl]). 

The properties described in Section 3.4 are valid for all variational integrators. In the following we present 
special properties of the variational integrators constructed in this work. 
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3.5 Time reversibility 

A further geometric property of Hamiltonian systems is the time reversibihty. It seems hkely to use numerical 
methods that produce a reversible numerical flow when they are applied to a reversible Hamiltonian system. 
Furthermore, the numerical solution has a long-time behavior similar to the exact solution (see [HLW02]), 
and thus time reversibility is a desirable property also for the variational integrators presented in this work. 
In the following, we repeat the definitions of time reversibility and the adjoint discrete Lagrangian (following 
[IILW02] and [MWOl]) und show, under which conditions the variational integrators (PsNrQu) are time- 
reversible. 

Definition 3.12 (Time reversibility and adjoint ([HLW02])). A numerical one-step method ^d is called 
symmetric or time-reversible if it satisfies 



d ) 



$Jj o $^ '' = id or equivalently <i>Jj = (<f> 
The adjoint method denoted by $^ is defined by 

A method is called self-adjoint if we have 

Thus, symmetric, time-reversible, and self-adjoint are equivalent notions. In the following, we show that the 
integrators {PsNrQ2rGau) and {PsNrQ2r — 2Lob) are self-adjoint and thus time-reversible methods. 

Definition 3.13 (Adjoint of the discrete Lagrangian ([MWOl])). The adjoint discrete Lagrangian L*^^ of the 
discrete Lagrangian is given by 

The discrete Lagrangian is self-adjoint if 

Ld{qlql,h)^Ll{qlql,h). 

The following well-known theorem (see e.g. [MWOl]) connects the adjoint of the discrete Lagrangian with 
the adjoint of the discrete Hamiltonian flow. 

Theorem 3.14. If the discrete Lagrangian Ld has a discrete Hamiltonian map F^^, then the discrete Hamil- 
tonian map of the adjoint discrete Lagrangian LJ equals the adjoint map, i.e. -Fl* = -^X ■ If the discrete 
Lagrangian is self-adjoint, then the method is self-adjoint. Conversely, if the method is self-adjoint, then the 
discrete Lagrangian is equivalent"^ to a self-adjoint discrete Lagrangian. 

A proof of this theorem can be found in [MWOl], Theorem 2.4.1. To show the time reversibility, we first 
show that the discrete Lagrangian is self-adjoint. 

Theorem 3.15. Let Ld be the discrete Lagrangian of (PsNrQu) with symmetric quadrature formula (bi, Ci)l^i 
and interpolation polynomial qd with symmetrically distributed control points (di)f=o, i.e., bi — br+i-i, 
Ci + Cr+i-i = 1, i = 1, . . . , r, and di = 1 — ds-i, i = 0, . . . , s. Then Ld is self-adjoint. 



Proof. We have 



qd{-t;qo,-h) = '^qolu,s{t/h) ^qd{t;qQ,h), (3.8) 

i/=0 

1 " 

qd{-t-qo,-h) = —-^qQl^^sit/h) ^ -qd{t;qa,h). (3.9) 



-h 



^ Two discrete Lagrangian are equivalent if their discrete Hamiltonian map are the same. 
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Due to the symmetry of the control points di it follows for the Lagrange-polynomials l^^^s with t e [0, 1] 

(1 - r) - d, yj (-r) + (1 - dO 



'-('-' - n ^^x^^ n j^ 



( -r) + 4_^ _ -j-r r - ds-i 



n v^,7:'\ - n 

With go = (go, ••■, Qo), Qo = Qi, and go = (<7o, •■-, 1o) we have for t e [0, /i] 



qd{t,qo,h) = ^ go ""/j,, J - j = ^ go^5-fc,s 



= E 9o^fc,. ( 1 - ^) = E ^o^fc.s ( ^ 

fc=0 ^ ^ fc=0 ^ 

= qd{h-t;qo,h) 
and by taking the time derivative we have 

qd{h-t;qo,h) = -qd{t;qo,h). 

Thus, together with (3.8)-(3.9) we obtain 

qd{-'t;qo,-h) = qd{t;qo,h) ^qd{h-t;qo,h), 
qd{-t;qo,-h) = -qd{t;qo,h) ^qd{h-t;qo,h). 

By substituting these expressions in the adjoint discrete Lagrangian we have with the symmetry of the 
quadrature formula (6i,c,i)r=o 

r 

-Ld{qi,qQ,-h) = ext - {-h)'^biL{qd{-c,h;qo,-h)),qd{-c,h;qo,-h)) 

i/e{i,...,s-i} 

r 

= ext h^biL{qd{h- Cih;qo,h)),qd{h- c^h;qo,h)) 

i'e{i,...,s-i} 

r 

= ext h'^br+i-iL{qd{cr+i^ih;qo,h)),qd{cr+i-ih;qo,h)) 

J/G{l,...,s-1} 

= ^d(go°,g?)- 



Theorem 3.16. T/ie discrete Hamiltonian maps of {PsNrQ2rGau) and {PsNrQ2r~2Lob) are self-adjoint 
and thus time-reversible. 

Proof. Since the Gauss and the Lobatto quadrature are symmetric and since the control points of qd are 
chosen symmetrically, it follows with Theorem 3.15 that the discrete Lagrangian of {PsNrQ2rGau) and 
{PsNrQ2r — 2Lob) are self-adjoint. The statement follows with Theorem 3.14. ■ 
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3.6 Linear stability analysis 

In the following, we investigate the stability properties of the constructed variational integrators. We restrict 
ourselves to a linear stability analysis. Following [LR04], we consider the Lagrangian of a harmonic oscillator 

(3.10) 



r / -N -'^ -2 1 2 2 

L{q,q) ^ ^q " 2^ * 




with uj,q eR. The Hamiltonian equations read 




p=-uj'^q, q^p 




and the exact solution is given by 




/ p{t) \ f cos (Ljt) -sm{u}t) \ / p{0) 
y Luq{t) J y sin (ujt) cos (wi) J \ tjq{0) 


]-M ( P(«) 


with det{M^) — 1. The eigenvalues of M,^ are Ai^2 = e"*"™* and thus, 


we have that |Ai_ 



: 1. 

By applying a variational integrator (PsNrQu) to the Lagrangian (3.10), the discrete Euler-Lagrange equa- 
tions form a linear system of equations such that the discrete Hamiltonian map F^^ of (PsNrQu) can be 
written as linear map 

Pi \ - M f Po 

with matrix M^^uj- 

Following [LR04] , we call a numerical solution asymptotically stable if the growth of the solution is asymp- 
totically bounded. A sufficient condition for asymptotic stability is that the eigenvalues of Mh^uj are in the 
unit desk of the complex plane and are simple if they lie on the unit circle. We investigate this property for 
selected variational integrators. 

1. For the midpoint rule {PlNlQ2Gau) we have 



^"''^^ ~ 4hu; h^u^-4 

Since Mh^^Mj^^ = Id, M^^^ is orthogonal with |A(Af^,^)| = 1. Thus, the midpoint rule is asymptoti- 
cally stable for all /ijW S K. 

2. The Stormer-Verlet method {PlN2Q2Lob) has the iteration matrix 

(1 _ h^ 

with eigenvalues 

Al,2 — 1 z — ± 




2 
thus, the method is stable for {hu)^ < 4 (cf. [LR04]). 

3. For the {P2N2QAGau) method we have 

M = I 'i''w* + 12/i2i^2_|_;^44 h^ui^ + l^h'^ui'^ + UA 

and the method is stable for all g, w G M since M/j^^j is orthogonal. 
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4. The {P2N3QALob) scheme results in the iteration matrix 



1,4, .4 oou2, .2 



Mh 



2h^Lj^+48 h^ui^+24 

hu{h''u,''-3eh^u;^+2as) /i4^4_^3^2^2^4g 

12h^ui^+288 2h'^ui^+48 



with eigenvalues 



Al.2 — 



h'^Lo^ - 22/i2a;2 + 48 ± hujVh^uj^ - 4.4. h^^^ + 576/i2cj2 _ 2304 



2/i2 w2 + 48 
The stability region is shown in Fig. 3.2. The integrator is stable for \hijj\ < 2\/2. 

|Al,2| 



-o- 



-o- 



hw 



Figure 3.2: Modulus of the eigenvalues of Mh^u: for the {P2N3Q4Lob) integrator. The method is stable for 

\huj\ < 2^2. 



5. The iteration matrix for {P3N3Q6Gau) reads 



M, 



h,uj 



/i°i.j°-264/i^t^-*+6480?i^a;^- 14400 

/l6w'i+24/l4a;4 + 720/t2(.c;2 + 14400 

2ihui[h'^uj^-10h'^u)'^+Qm) 
^ h<^u)<^+24h'^u}'^+720h^u}^ + lH00 



24/iw(/i*Ld*'-70?i^w^+600) 

/i6w6+24/t4(^4_|_720/i2cj2 + 14400 

/t'^Lj'^-264fe'^cj''+6480fe'^aj''- 14400 
h^u)f^+2Ah^Lj'^+720h'^u)'^ + 14400 



The scheme is again stable for all (7,0; G M due to the orthogonality of Mh^u- 
6. The integrator (PSNAQGLob) gives 



Mh.u. = 



~46h*w'"+840/i^tj^-1800 
h^u'i'+GOh'^ui'^ + 1800 



6?iw(/i* 



-40ft^ 



h300) 



h^ui'^+GOh'^ui^ + 1800 



hu;{h<^u,<'-144h''u;*+5760h^uj^-43200) h<' u^'' -92h^ lo" + lQ80h^ u^ ~3Q00 



24(/!,4t^4_|_60/i2^2_).1800) 



2/i4lj4 + 120/i2-i«2+3600 



The eigenvalues are given as 



Ai 



3600 - 1680x2 + 92a;'' 



x^ ± xn/xIO - 184x8 + 11820a;6 - 316800a;4 + 3456000a;2 - 12960000 



2x4 + 120x2 + 3600 
with X :— hu. In Figs 3.3 and 3.4 (zoom of Fig. 3.3) the stability region is shown. 



An interesting observation is that for the integrators {PsNsQ2sGau), s — {1,2,3} the iteration matrices 
-M^h,tj are orthogonal independent of the step size h and thus asymptotically stable. If this also holds for 
general integrators of this type has to be investigated in the future. 
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3- 




-7-6-5-4-3-2-101234567 



hw 



Figure 3.3: Modulus of the eigenvalues of Mh.uj for the {P3N4Q6Lob) integrator. 

|Ai,2| 4 
1.02 

1.01 

1.00 

0.99 

0.98 

3.10 3.11 3.12 3.13 3.14 3.15 3.16 3.17 3.18 

hw 

Figure 3.4: Modulus of the eigenvalues oi M^^^ for the {PiNAQQLoh) integrator (zoom). 




4 Variational error analysis 



In this section we analyze the order of accuracy of the constructed variational integrators (PsNrQu) for 
s,r,u G N. To this end, we introduce the concept of the variational order of a discrete Lagrangian (see 
[MWOl]). The main idea is the following: Rather than considering how closely the trajectory of the discrete 
Hamiltonian map matches the exact trajectory, we consider how closely the discrete Lagrangian matches 
the action integral. It was shown in [MWOl] and [PC09] that both order concepts are equivalent (see 
Theorem 4.2). 

To define the variational order of a discrete Lagrangian, we consider the exact discrete Lagrangian (see 
[MWOl]) that is defined to be 

Lf{qlqlh)= [\{q{t),q{t))dt 
Jo 

for sufficiently small h and close configurations qo and qi. q{t) is the unique solution of the Euler- Lagrange 
equations for L that satisfies the boundary conditions g(0) — (7° and q{h) — q^. 

Definition 4.1 (Local variational order ([MWOl])). A given discrete Lagrangian Ld is of order r if there 
exist an open subset Uy C TQ with compact closure and constant Cy > and hy > such that 

\\LMO),q{h), h) - Lf (<7(0), q{h), h)\\ < Cyh-+' (4.1) 

for all solutions q{t) of the Euler- Lagrange equations with initial condition {q,q) € Uy and for all h< hy. 
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An extended version of the following theorem can be found in [MWOl]. 

Theorem 4.2. Given a regular Lagrangian L and corresponding Hamiltonian H , then a discrete Lagrangian 
Ld is equivalent'^ to a discrete Lagrangian of order r if and only if the discrete Hamiltonian map for Ld is 
of order r. 

In the following section we want to determine the error defined by (4.1). To this end, we assume the solution 
q{t) to be smooth enough, in particular we assume q € C*"'"^'''+^([0,T], Q)"^ and define by 



(t; (qK/i)):^o), /i) = ^ gK/j)Z,, J| 



Id 

the polynomial qd of degree s on [0, h] that interpolates a solution of the Euler-Lagrange equations q{t) in 
t = dih, i = 0, . . . , s. By the interpolation error for the Lagrange interpolation, we have that 

q{t) = qd{t) + 0{h'+^), (4.2a) 

q{t) = qd{t) + 0{h') for all t e [0, h]. (4.2b) 

In particular, at the control points, we obtain 

qd{dih;{q{d^h))l^f^),h) = q{d,h), (4.3a) 

qd{d^h■{q{d,h))Uo),h) = q{d^h)+0{h') for i - 0, . . . , s. (4.3b) 

For the variational error analysis of the integrators (PsNrQu) we have to consider the discrete Lagrangian 

r 

Ld{q{0),q{h),h) ^ h^ hL{qd{cih; {q{d^h))l^Q, h),qd{cth; {q{d^h))l^Q, h)), 

where q(t) satisfies the Euler-Lagrange equations. The interpolation error (4.2) allows us to obtain the 
following theorem of the minimal order of accuracy of the discrete Lagrangian. 

Theorem 4.3 (Minimal order of (PsNrQu)). Let L be a Lipschitz-continuous Lagrangian in both variables 
and q{t) the unique solution of the Euler-Lagrange equations with q(Q) — qg and q{h) — qj*. Let the trajectory 
q(t) he approximated by a polynomial qd{t) of degree s and let the action integral be piecewise approximated 
by a quadrature formula {bi,Ci)l^i of order u. Then, the corresponding discrete Lagrangian is at least of 
order min(s, u). 

Proof. A comparison of the discrete and the exact discrete Lagrangian leads to 

r 

Ld{qiO),q{h)) = h y^ biL{qd{cih),qd{cih)) 

1=1 

= h\y^hMq{c,h),q{c^h))+0{h')\ 

r 

= hY,b^L{q{c,h),q{c^h))+0{h'+^) 

i=l 

h 
L{q{t),q{t))dt + 0(;i"+l) + 0{h'+^) 







= Lf (9(0), g(/l)) + 0(/i'nin("+l,-+l)) 



^Two discrete Lagrangian L^ and L^ are equivalent if they have the same discrete Hamiltonian map, i.e. F]^^ = Fj^ 
^C*''°([0,T], Q) is the space of functions q : [0,T] — > Q virhich are k times continuously differentiable on (0,T) and whose 



fc-th derivative is Lipschitz continuous on [0,T]. 
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A proof of this theorem can also be found in [LSlla]. However, in contrast to [LSlla], we emphasize 
that min(s,w) is only the minimal order of the discrete Lagrangian and the associated variational integra- 
tor, respectively. In the following we want to determine the maximal order of the variational integrator 
(PsNrQu). In Section 3.3 we assumed that the discrete Lagrangian is regular to obtain a well-defined dis- 
crete Lagrangian flow. In [IIL12] it is shown that the well-posedness depends on the order of the quadrature 
rule that is used to approximate the action integral. In particular, it is shown that for a Lagrangian of the 
form L = ^q^Mq—V{q) with M symmetric positive-definite and VV^ Lipschitz continuous a unique solution 
of the internal stage equations (3.4) exists if the used quadrature rule is of order at least 2s — 1 (with s the 
degree of the interpolation polynomial). Thus, to ensure a well-defined discrete Lagrangian flow, we make 
the following 

Assumption 4.4. For the variational integrators (PsNrQu) we have u> 2s — 1. 

For the Gauss and Lobatto quadrature this means that we have to choose r > s and r > s -|- 1, respectively, 
since this yields quadrature rules of order at least 2s. Note that the variational integrator {PsNsQ2s — 
2Lob) which yields the standard Lobatto IIIa-IIIB partitioned Runge-Kutta method does not satisfy the 
assumption, however, the condition m > 2s — 1 is only a sufficient not a necessary condition for uniqueness 
of solutions. 

For the error analysis we assume a regular C^'^ Lagrangian L : TQ -^ M with fc > 2s + 1 and n-dimensional 
configuration manifold Q and consider the Taylor expansion of Lci{q{0),q{h)) in {q,q) = (g(0),g(0)) which 
reads 

LMO),QW,h) 

r 

= h^ b^L{qd{c^h; {q{dt,h))l^g), h), qd{c^h; {q{d^h))l^g), h)) 

i=l 
- '^tb.tl±('Ag^Uqdic.h)-qyAUc.h)-qr^)+Oih^^^^^^ 

For brevity, in the third line we just write qd{cih) for qd{cih; {q{di,h))f,^Q), h). The multivariate derivatives 
of the Lagrangian L are given by the multilinear functions 

j times k—j times 

ioT j — 0, . . . ,k and fc = 0, . . . , 2s. 

Theorem 4.5. Let L be a C^'^ Lagrangian (k > 2s + 1) and let qd{t) be a polynomial of degree s that 
interpolates q{t) in djjh, v — Q^ . . . ,s, with control points {d^)l^Q given by the quadrature points of the Lobatto 
quadrature (w^, (i,y)f,^Q w.r.t. the interval [0,1]. Let (q^q) — (q(0),g(0)) and let {bi,Ci)l^i be quadrature 
formula w.r.t. the interval [0, 1] of order u satisfying Assumption 4.4- Then it holds for all k = 0, ... ,2s 



'^E ^'C^M((«'^(^^^) - 9F' iAdc^h) qf-^) 



' dqWq''~i 



= hjlb.^t^JP^Uc^h) - q)\ {q{c,h) - q)''-^) + 0(/,-in(2^+i.«+i) 



=1 



dqidq^ ^ 



Proof. For brevity we assume {q{Q),q{Q)) — (0,0). This is no restriction because all steps of the proof can 
be performed with the polynomials qd{t) = Q_d{i) — 9(0) and (?d(i) = 'Zd(i) — 9(0) instead of qd{t) and gd(i), 
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respectively. We compute 



hjlb.^^^{{qMh)y,{Uc^h)f~') 



i=l 
h Qfc 






. dqWq^ 



d''L{q,q) 



= h'^uji ^ — ^ — ■ qd,n{dih)- ■■qdAj{dih) ■ qdjiidih)- ■■qd,h-A'^i^) 

I /n/Lmin (2s+l,M+l)\ 
h---lk-j — ^ 

[qii {dth) + [qd^i^ {dih) - qi^ (dih))) ■ ■ ■ {qi^_. (d^h) + [qd.h-j (dih) - qi^^^ (dih))) 

I ^/Lmin (2s+l,u+l}\ 

'^ " d^'Lia a) 
= /i ^ Wi ^ g ...g g. .. a- 9*1 (fii/i) • • ■ ft, (rfz/i) • 9(1 {d^h) ■ ■ ■ qi^_^ (d^h) 

d^Liq q) ^~^ '^~"' 

+ ^^ ^ ^^^ ^-^ Qiiidih) ■■ -qi^dih) ■ '^{qd,ia,{dih) - qi^{dih)) J]^ qift{dih) 

(Jqii' ■ ■ clqijiOqi-^^- ■ ■ oqi^_. ^^^ ^^^ 

I ^/Jjinin (2s+l,u+l)\ 

= /iVwi V — ^ — q,^{dih) ■ ■ ■ q^{d,h) ■ qi^{d,h) ■ ■ ■ qi^_{d,h) 

U .1.^=1 %i-- ■%,,%! ■■•%.-, 

_|_/n/'J,min (2s+l,u+l)-i 
afc , 



= /»E^'|^M((9(^'Mr,(9(rf.M)'-o+0(/^"""^'^+''"+'^) 



'i=0 



dqWq^ 






j=l 



dqidq^ 



where the first and the ninth equahties follow from the order u of the quadrature formula {hi^CiY^^i. The 
second and the eighth equalities follow since the Lobatto quadrature (w^,di,)j^^Q is of order 2s (cf. [Hil87]) 
The fourth and fifth equalities result from the equations (4.3) and by summarizing terms of higher order. 
The sixth equality follows from the following calculation: As last term in the sum {a — k — j) we obtain 
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h 



'dqi, ■■■dqi.,dqi, ■ ■ ■ dqi^_ 
d''L{q,q) 



dq^^---dqi^,dqi^---dqi^_^ 
d'^L{q,q) 



Qh (dih) ■ ■ ■ q,. (dih) ■ qi^ {d,h) ■ ■ ■ qi^_._-, (dih) ■ {qd.ik-j (dih) - qi^_^ (d^h)) 
q^, (*)••• q^, (t) ■ qi, (i) • • • g,,_^-_, (t) • ((?d,,,_, (t) - qi,_^ {t)) dt + 0{h^'+') 



%i ■■■dqi^,dqi^ ■■■dqi^_^ 
d''L{q,q) 



Qti (*)■■• qt, (t) ■ qh (*)■■• qi^-,-i (t) ■ {qdM-, W - qh-, (*)) 



(4.4) 



d 
di 



dqii ■■■dqt.,dqi^ ■■■dqi^_. 



qiAt) ■ ■ ■ qi,{t) ■qh(t)---qh-j-i{t) 



{qdj,^At)-i'.-At))dt + o{h'^+') 



^hj:. 



dt 



d''L{q,q) 



%i ■■■dq^^,dqi^ ■■■dqi^_^ 



q.iAdih)- --qi-idih) ■ qi-,{dih) ■■ ■ qi^^-_^{d,h) 



■{qd^,^Ad^h)-ql,_Ad^h))+0{h''+') 

where the the first and third equalities follow from the order of the quadrature formula {uj,^,d,y)l^Q and 
qd{,Q) = 9(0), qd{h) = q{h), the second equality by using integration by parts and the fourth equality with 
equation (4.3a). For the remaining terms {a — 1, . . . , k—j — 1) the same computations can be performed. ■ 

We finally state the following 

Theorem 4.6 (Order of (PsNrQu)). Let Ld the discrete Lagrangian with quadrature formula {hi,Ci)\^i 
of order u satisfying Assumption 4-4 ^.^d qd{t) an interpolation polynomial of degree s with control points 
{diy)f,^Q given by the Lohatto quadrature {uJi^,d,y)l^Q. Then L^ is of order min {2s, u). 



Proof. By Taylor expansion oi Ld{q{Q),q{h)) in {q,q) = (q(0),g(0)) and with Theorem 4.5 it follows 
Ld{q{0),qih),h) 

r 

= h^biL{qd{cih; {q{d^h)yi^Q),h),qd{cih; {q{d^h)yi^„),h)) 



k \ d^L{q,q) 



■r 2s k 

^ ^ ^' ^ fc! ^ V J J dqldq'^-l 
j=i fc=o i=o ^ ■" ^ ^ ^ 

+ 0{h^'+^) 

r 2s ^ k 



{qd{cih) - qy,{qdicih) - q) 



k-j 



^E^^E^E 



i=i fc=o j=a ^ ■' ^ ^ ^ 



k \ d'^Liq^q) 



(qic^h) ~q)\{q{crh) - q) 



k-3 



I ^ / i^min(2s + l .u+l) \ 

r 
/l^&,L(q(c,/l),g(c,/l))+0(/l--(2s+l,n+l)) 



rh 



f L{q{t),q{t))dt + C)(/imi„(2s+l,„+l))^ 
Jo 



Remark 4.7. One crucial assumption in the proof of Theorem 4-5 is that the control points {d^)'l^Q of 
the interpolation polynomial are given by the quadrature points of the Lobatto quadrature. In our numerical 
examples we observed that also a different choice of the inner control points {d,y)1'Z-^ leads to the same 
convergence rate. This might be due to the presumption that the extremizer of the discrete Lagrangian 
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constructed with a polynomial of degree s with control points = da < di < ■ ■ ■ < ds~i < ds = 1 and the 
extremizer of the discrete Lagrangian constructed with a polynomial with control points = do < di < ■ ■ ■ < 
ds-i < ds = 1 give the same polynomial qd(t) as solution but formulated based on different control points. 
This has to be investigated in the future. 

Theorem 4.8 (Maximal order of (PsNrQu)). The maximal order of the variational integrator (PsNrQu) 
is 2s. 

Proof. We assume that the Lagrangian {PsNrQu) is of order 2s+l and let (6^, Ci)l^i be the used quadrature 
rule with arbitrary order of at least u = 2s + l. Let q{t), t e [0, h], be the exact solution of the Euler-Lagrange 
equations and qd{t), t G [0, h], the polynomial of order s that interpolates q{t) in i e {doh, . . . ,dsh} such 
that 

qdid,h)^q{0) = q{d,h)-q{0), for alH = 0, . . . , s. (4.5) 

The control points (di)f^o °^ Qdit) can be chosen freely except for rfp = and dg = 1. Further, let (a^, di)f^Q 
be a quadrature formula with the control points di of qd{t) as quadrature points. Its degree of accuracy is 
of at least s, since it uses s + 1 quadrature points (cf. e.g. [Hil87, SK06]). 

Since the Lagrangian (PsNrQu) is of order 2s + 1, the coefRcient comparison in Theorem 4.6 shows in 
particular for A: = j = 1 that the following has to hold 

^E ^»l^('?(0)' '?(0))(9'^(^^^) - 9(0)) = hJ2h~{q{fS).,m){q{c,h) - q{0)) + 0{h^^+^). (4.6) 

Since ^{q{0),q{0)){qd{t) — q{0)) is a polynomial of degree s, it is integrated exactly by the quadrature rules 
{bi,Ci)l^i and {ai,di)f^Q and we have 

'^81 ^ F)T 

^{q{0),m){q{t) - q{0)) dt = /» J] 6i^(g(0), 9(0))(g(c,/z) - q{0)) + 0{h'^+^) 



=1 



dq 



^'= hY,b~{q{Q)Am{qd{c.h) ~ q{0)) + 0{h'^+^) 

, ^{q{0),q{0)) {qd{t) ~ qm dt + 0{h'^+') 
Oq 

= hY,a~{q{0),m)iqd{d^h)-q{0))+O{h^^+') 
1=0 9 

^ = ^ /i^a,|^(g(0),9(0))(g(d,;i)-9(0)) + O(/^2^+2)^ 

It follows that the quadrature formula {ai,di)f^Q integrates the function '^{q{0),q{0)){q{t) — q{0)) exactly 
up to an error of ©(/i^^^^), and thus, is of order 2s + 1. This is in contradiction to dp = and dg = 1 because 
by fixing these control points the order of the quadrature formula is at most 2s (cf. e.g. [Hil87, SK06]). Thus, 
the maximal order is 2s which is indeed possible by Theorem 4.6. ■ 

By Theorem 4.6 and 4.8 we can conclude that for a variational integrator based on an approximation space of 
degree s polynomials the convergence order 2s is possible. This maximal order is reduced if the quadrature 
rules used for the approximation of the action is not accurate enough. Thus, to guarantee the maximal 
convergence order we have to choose u > 2s. 

Remark 4.9 (Discussion on efficiency). Theorem 4-6 indicates which combination of quadrature rule and 
polynomial degree leads to lowest computational effort for a given order. As will be demonstrated in Section 5, 
the computational effort grows for an increasing number of quadrature points and control points. A polynomial 
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Table 2: Variational order of {PsNrQ2rGau) 
s=4 s=5 s=6 



r\s 


s = 2 


s = 3 


r = 2 


4 




r = 3 


4 


6 


r = 4 


4 


6 


r==5 


4 


6 


r = 6 


4 


6 



8 10 

8 10 12 

of degree s has s + 1 control points which leads to s + 1 Euler- Lagrange equations that have to be solved for. 
Since the computational effort grows if more Euler- Lagrange equations have to he solved for and since the 
order of the integrator is niin(2s, u), a reasonable choice for the polynomial degree is half of the order 
of the quadrature formula, i.e., u — 2s. This guarantees a minimal number of control points without an 
order reduction of the variational integrator. For the Gauss quadrature rule we thus choose the integrator 
(PsNrQ2rGau) with r = s. As discussed in [MWOl, HLW02], the Lobatto quadrature is advantageous for 
stiff systems and Theorem 4-6 indicates the choice r — s + 1 for the variational integrator {PsNrQ2r — 2Lob) 
to guarantee an order 2s. 

Remark 4.10 (Discussion on (PsNsQu) ([MWOl])). The integrator {P sN sQ2sGau) corresponds to the 
collocation Gauss-Legendre rule which is known to be of order 2s. Theorem, 4-6 together with Theorem 4.2 
provide the same result. The integrator {PsNsQ2s — 2Lob) is the standard Lobatto IIIA-LIIB partitioned 
Runge-Kutta method and of order 2s — 2 (as known from theory on Runge-Kutta integrators). Although 
this integrator does not satisfy Assumption 4-4; ihe statement of Theorem 4-6 is still true in this case. In 
particular. Theorem 4.6 shows that the order of {PsN sQ2s — 2Loh) does not decrease if a polynomial of degree 
s — 1 instead of s is used, or, the order increases to 2s if a Lobatto quadrature with s + 1 quadrature points 
is used. In general, it is only efficient to choose the number of quadrature points equal to the polynomial 
degree if the Gauss quadrature is used because only in this case the order of the integrator (PsNsQu) is not 
reduced by the quadrature formula. If no Gauss quadrature rule is used, we have a quadrature order u < 2s 
and with Theorem 4-6 we obtain a convergence order of u. A polynomial of degree [u/2] < s leads to the 
same convergence order, however to s — [m/2] less discrete Euler- Lagrange equations. 



5 Numerical examples 

In this section the presented variational integrators are investigated numerically and the theoretical results 
on the order of accuracy as stated in Theorem 4.6 as well as the preservation properties are evaluated. To 
this end, we consider two examples, the harmonic oscillator and the Kepler problem. We restrict to the most 
promising integrators w.r.t. high accuracy and low computational effort as discussed in Remark 4.9, i.e., the 
polynomial degree is smaller or equal to r and the quadrature rule in use is of Gauss or Lobatto type. By 
Theorem 4.6 we know that the variational integrators {PsNrQ2rGau) and {PsNrQ2r — 2Lob) are of order 
min (2s, 2r) and min (2s, 2r — 2), respectively. The orders for the different versions of these integrators are 
listed in Table 2 and 3, respectively. The discussion in Remark 4.9 as well as the following numerical results 
show that the relations r = s and r = s + 1 are the best choices for a Gauss and a Lobatto quadrature, 
respectively. The implementation is performed as described in Remark 3.5. 



5.1 Harmonic oscillator 



Consider the two-dimensional harmonic oscillator with mass being equal to one and the Lagrangian L(g, q) 
\q^q — \q^ q with (7,(7 G M^. The Euler-Lagrange equations are 

m = ~q{t). 
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Table 3: Variational order of {PsNrQ2r - 2Lob) 

■8 = 4: s = 5 s = 6 
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r = 4 
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r = 5 
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r = 6 
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6 



10 



10 



By the Legendre transform we obtain the Hamiltonian system 

m^P{t), p{t) = -q{t). 



(5.1) 



The total energy of the system is given by the Hamiltonian H{q,p) = ^p p + ^q q with q.p e M . Let 
{q'^,p'^){t) denote the exact solution consisting of configuration and momentum of the Hamiltonian sys- 
tem (5.1). For the error calculations we use the global error determined by 







max qi. 
ke{o,...,N} ' " 
ie{i,2} 



qKhk)l 



max \p1 . 

ke{Q,...,N} 

ie{i.2} 



■pKhk)\ 



(5.2) 



with step size h and the index i denotes the components of the configuration and momenta, respectively. 
ilk'Pt)k=o i^ ^^^ discrete solution computed by a variational integrator with qJv — 9^-1 ^-^d p'^ = Pn-i- ^^ 
Figs 5.1 and 5.3 the global error for the different variational integrators (PsNrQu) is shown in dependence 
of the step size h. The numerically determined order is given in Figs 5.2 and 5.4. A comparison with Tables 2 
and 3 shows that the orders match with the analytically determined values from Theorem 4.6. 
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Figure 5.1: Harmonic oscillator: Log-log plot of the error for the configurations q and the momenta p 
(superimposed) by step size h; with the use of variational integrators [P sNrQ2rGau)] divided into four 
subplots, separated by applied Gauss quadrature. 

In Fig. 5.13 the run-time compared with the global error is given for the different integrators {PsNrQ2r — 
2Lob) with s = {1,...,6} and r = {3, ...,6} and for different step sizes 
h e {1,0.5,0.25,0.125,0.1,0.0625,0.03125}. The step size is reduced until the error is below 10"^°. The 
integrators {P2N3QALob), {P2NAQ6Lob), {P2N5Q8Lob), and {P2N6Q10Lob), which are of order four, 
demonstrate that for an increasing number of nodes also the run-time increases as explained in Remark 4.9. 
The same behavior is also observable for the other integrators. Further, we explained in Remark 4.9 that a 
higher polynomial degree leads to an increasing number of Euler-Lagrange equations. In Fig. 5.13 it is shown 
that this leads to an increasing run-time (compare for example {P2N'iQALob) and {PZNiQALob) which are 
both of order four; or {PbNQQlQLob) and (P6A^6Q10Lo5), which are of order ten). 
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Figure 5.2: Harmonic oscillator: Order of the variational integrators {P sN rQ2rGau) with respect to the 
configurations q (left bar) and the momenta p (right bar). The title gives the order of the applied Gauss 
quadrature which is 2r. The x-axis indicates the degree of the used polynomial. 



Lobatto quadrature of order 4 Lobatto quadrature of order 6 Lobatto quadrature of order 8 Lobatto quadrature of order 1 
10= I . 1 ^fy' I . 1 10= I : 1 10* I ' 1 




Figure 5.3: Harmonic oscillator: Log-log plot of the error for the configurations q and the momenta p 
(superimposed) by step size h; with the use of variational integrators {PsNrQ2r — 2Lob); divided into four 
subplots, separated by applied Lobatto quadrature. Note that {Ps— lNsQ2s — 2Lob) and {PsNsQ2s — 2Lob) 
are of the same order. 

Hence, Fig. 5.13 underpins Remark 4.9 and justifies the conclusions made there. That is, the integrators 
{PsNsQ2sGau) and (PsNs + lQ2sLob) are the most efficient integrators with view of run-time per order. 

For a clearer illustration we omit in Fig. 5.14 the less efficient integrators which are displayed in Fig. 5.13 
and include the most efficient integrators constructed with the Gauss quadrature. It can be read off which 
of these integrators provide the desired accuracy with shortest run-time. Larger s means longer processing 
time but also higher order. Notice that {PlNlQ2Gau) and {PlN2Q2Lob) are the midpoint rule and the 
Stormer-Verlet method, respectively. 



5.1.1 Conservation of angular momentum and long-time energy behavior 



Positioning the system of the two-dimensional harmonic oscillator in the {x, y)-plane of the three-dimensional 
space, the z-component of the angular momentum 

i{v,(i) = -piq2+P2qi 

is a conserved quantity. This follows from Noether's theorem because the Lagrangian of the two-dimensional 



harmonic oscillator L = ^q q— i^q q is invariant under the group of rotations 5*0(2) = {B e 



\B'^B 
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Figure 5.4: Harmonic oscillator: Order of the variational integrators {PsNrQ2r — 2Lob) with respect to the 
configurations q (left bar) and the momenta p (right bar). The title gives the order of the applied Lobatto 
quadrature which is 2r — 2. The x-axis indicates the degree of the used polynomial. 

Id, det(i3) = 1}. With the linearity of i?„ G 50(2) and Remark 3.11 it follows that all variational integrators 
(PsNrQu) conserve the ^-component of the angular momentum. In Fig. 5.5 it is shown that the z-component 
of the angular momentum, if the variational integrators (P2iV3Q4Lo6), {PiNAQQLoh) and (PANbQSLob) 
are used, is preserved up to an error less than 10""'^'*.'* In Fig. 5.5 (right) the behavior of a Runge-Kutta 
method of order four is included which is not symplectic nor momentum-preserving. Thus, Noether's theorem 
does not apply and the angular momentum is not preserved. In Section 3.4 we introduced the good long-time 
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Figure 5.5: Harmonic oscillator: Left: Error of the z-component of the angular momentum for the variational 
integrators (P2iV3Q4io6) (dots), {P3NAQ6Lob) (crosses) and (PiNbQSLob) (squares) with step size h = 
0.5. Right: Same plot as on the left with non-variational integrator included (Runge-Kutta method of order 
four with step size h = 2^^ (red solid)) . 

energy behavior of symplectic integrators, in particular of variational integrators. In Fig. 5.6 the error of 
the total energy of the harmonic oscillator simulated by different integrators is shown. While the use of the 
variational integrators {P2>N AQQLob) and (PANbQSLob) leads to an oscillating but stable energy behavior, 
the use of a nonsymplectic Runge-Kutta method of order four clearly shows an energy drift. 



*Note that the accuracy is limited to machine precision and the accuracy of the applied Newton method to solve the discrete 
Euler-Lagrange equations. 
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Figure 5.6: Harmonic oscillator: Error of the energy for the integrators {PiN AQQLoh) with step size h = 0.25 
(crosses), (PANbQSLob) with step size h = 0.5 (dots) and for a Runge-Kutta method of order four with step 
size h = 2~^ (solid line). 

5.2 Kepler problem 

The two-body problem, also known as Kepler problem, is to determine the motion of two point particles with 
masses mi, 7712 € K. By assuming mi = 1 with gravitational constant 7 and k — jmim2 G M we construct 
the Lagrangian of the Kepler problem 



and the Hamiltonian 



L{qA) 



Hiq,p) = 



1 .T. 
(7 (7 + 


k 




2? 9 + 
1 J. 


k 


9 
92 


2^ ^ 


V<i! + 


4 



(5.3) 



with q = {qi,q2)^ ,q,p G M^ . The Hamiltonian equations provide the system of differential equations which we 
want to solve with respect to the initial condition {qo,po). For the simulations we set k — 1.016895192894334, 
go ~ (5, 0)-^ and po = (0, 17)-^ since this results in a motion where the mass mi describes an elliptic orbit 
around mass 7712 with period T = 5.0. Therefore, after ^ steps of step size h the simulated mass should 
return the fifth time to the initial point {qo,po). 

As error we compute the maximal difference between the given initial value and the value after N — "^ steps 
of integration given as 

and max \p^j, 



niax \q^ 
*e{i,2} 



qa,: 



»e{i,2} 



^%, 



Pa,: 



for configuration and momentum, respectively. The error is shown in Fig. 5.7 and Fig. 5.9 for the different 
variational integrators {PsNrQu) and for step sizes ft. G {1, 0.5, 0.25, 0.125, 0.1, 0.0625, 0.03125}. The numer- 
ically determined order is given in Fig. 5.8 and Fig. 5.10 and coincides nicely with the analytically predicted 
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order from Theorem 4.6 as given in Tables 2 and 3. It can be observed in the plots that the error of the 
integrators decreases not far below 10~^°. Since the iteration is implicit and has to be solved by a Newton 
method, the accuracy is limited by the machine precision and the used solver (we used f solve implemented 
in Matlab). 
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Figure 5.7: Kepler problem: Log-log plot of the error for the configurations q (crosses) and momenta p (dots) 
with step size h for the integrators {PsNrQ2rGau)] divided into four subplots, separated by applied Gauss 
quadrature. In the third subplot Ps denotes (PsNAQSGau) and in the fourth Ps denotes {PsN5Q10Gau). 
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Figure 5.8: Kepler problem: Order of the variational integrators {PsNrQ2rGau) with respect to the config- 
urations q (left bar) and the momenta p (right bar) . The title gives the order of the applied Gauss quadrature 
which is 2r. The a:-axis indicates the degree of the used polynomial. Pb gives the numerically determined 
order with use of the first six values whereas P5* gives the numerically determined order with use of the 
first five values. 



5.2.1 Conservation of angular momentum and long-time energy behavior 

Since the Lagrangian (5.3) of the Kepler problem is invariant under the group of rotations 5*0(2) = {B £ 
]j2x2 I qTq ^ jjj^ det(B) = 1}, the angular momentum 



np,q) 



-piq2 +P2qi 



is a conserved quantity in the system. The angular momentum is (up to numerical accuracy) also preserved 
in the discrete solution using the variational integrators (PsNrQu) (cf. Remark 3.11) as shown in Fig. 5.11 
for the integrators (P2iV2Q4Gau), {P3N3Q6Gau) and (PiNAQSGau). However, using a Runge-Kutta 



27 



Lobatto quadrature of order 4 Lobatto quadrature of order 6 Lobatto quadrature of order 8 Lobatto quadrature of order 1 




10" 



10 



P1N5Q8 



P2N5Q8 



P3N5QI 




10" 



10 





: 


PI -, 


4j-^-:M 


P2 ■/ 


7i 


PsJ^ 


// 


/p4j/ 


/P5 


■'r' 


P6 



10"' 



10 
h 



10" 



10"' 



10" 
h 



10" 



Figure 5.9: Kepler problem: Log-log plot of the error for the configurations q (crosses) and momenta p (dots) 
with step size h for the integrators {PsNrQ2r — 2Lob); divided into four subplots, separated by applied 
Lobatto quadrature. The values of {Ps — lNsQ2s — 2Lob) and {PsNsQ2s — 2Lob) are superimposed, they 
are of the same order. In the last subplot Ps denotes {PsN6Q10Lob). 
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Figure 5.10: Kepler problem: Order of the variational integrators(PsA^rQ2r — 2Lob) with respect to the 
configurations q (left bar) and the momenta p (right bar). The title gives the order of the applied Lobatto 
quadrature which is 2r — 2. The a;-axis indicates the degree of the used polynomial. 

integrator of order four, angular momentum is not preserved anymore as shown in Fig. 5.11 (right). The 
error in the energy is given in Fig. 5.12. As for the harmonic oscillator, due to the symplecticity of the 
variational integrators the error of the integrators {P3N3Gau) and (PANiGau) is oscillating but bounded 
whereas the Runge-Kutta solution exhibits an energy drift. 



6 Conclusion 



In this work variational integrators of higher order have been constructed by following the approach of 
Galerkin variational integrators introduced in [MWOl]. Thereby, the solution of the Euler-Lagrange equa- 
tions is approximated by a polynomial of degree s and the action by a quadrature formula based on r 
quadrature points. The restriction to r = s quadrature points (as assumed in [MWOl]), which leads to 
symplectic partitioned Runge-Kutta methods, is dropped. For the resulting methods the order of accuracy 
has been determined numerically as well as analytically. It has been shown that the order of accuracy can 
be increased by adapting the number of quadrature points to the polynomial degree. In particular, if the 
Lobatto quadrature is used, the choice of r = s -|- 1 leads to an integrator of order 2s which is of two orders 
higher compared to an integrator with r = s quadrature points (order 2s — 2). Thus, the ideal ratio between 
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Figure 5.11: Kepler problem: Left: Error of the angular momentum for the variational integrators 
{P2N2Gau) (dots), {PiNiGau) (stars) and (P47V4Gau) (squares) with step size h = 0,25. Right: Same 
plot as on the left with non-variational integrator added (Runge-Kutta method of order four with step size 
h = 2-6 (red solid)). 
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Figure 5.12: Kepler problem: Error of the energy for {P3N3Q6Gau) with step size h — 0.125 (squares), 
(PANAQSGau) with step size h ~ 0.25 (diamonds) and for a Runge-Kutta method of order four with step 
size h — 2~^ (solid line). 

the number of quadrature points r and the polynomial degree s could be determined for different quadrature 
rules and it has been shown that 2s is the maximal possible order of the constructed variational integra- 
tors. The analytically determined orders as well as the structure-preserving properties such as symplecticity, 
momentum-preservation and good long-time energy behavior have been demonstrated by numerical exam- 
ples. In addition, for symmetrically distributed control points of the polynomial, the variational integrators 
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{PsNrQ2rGau) and {PsNrQ2r ~ 2Lob) have been shown to be time-reversible. Furthermore, a hnear sta- 
bility analysis has been performed for selected integrators and stability regions have been determined. It has 
been shown that the integrators {PsNsQ2sGau) for s — {1,2,3} are asymptotically stable independent of 
the step size h. If this also holds for general integrators of this type has to be investigated in the future. In 
addition, it would be of interest to relax the condition on the control points in Theorem 4.5 (cf. Remark 4.7) 
since in the numerical experiments we observe the same order also if the control points do not correspond 
the the quadrature points of the Lobatto quadrature. 

In the future, the order of accuracy by using different interpolation methods, e.g. Hermite polynomials 
as introduced in [LSllb] to ensure higher smoothness at the time nodes, has to be investigated. The 
application of the constructed higher order variational integrators to holonomic and nonholonomic integrators 
is straightforward, however, a careful error analysis has to be carried out to see if the predicted orders also 
hold for these systems. Currently, the approach is used for the optimal control of mechanical systems 
([CJ0B12]) and numerical results confirm that also the adjoint resulting from the necessary optimality 
condition inherits the same order of the variational scheme as it also does for symplectic partitioned Runge- 
Kutta methods (cf. [OBJMll]). Another topic of interest is the construction of time-adaptive variational 
integrators since naive time-adapting strategies destroy the structure-preserving properties (see e.g. [LR04]). 
The higher order integrators could be applied to adapt the order rather than to adapt the step size, e.g. a 
higher order scheme can be deployed if higher accuracy requirements have to be matched. Furthermore, for 
systems involving slow and fast time scales (see e. g. [L0B13]) variational integrators of different orders for 
the different subsystems can be used to increase the efficiency of the simulations. Thereby, an investigation 
regarding preserved quantities and long-time behavior is essential. 
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Figure 5.13: Harmonic oscillator: Global error against run-time for the integrators {PsNrQ2r — 2Lob) and 
different step sizes h. 
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Figure 5.14: Harmonic oscillator: Global error against run-time for the integrators {PsN sQ2sGau) (dotted 
line) and {PsNs + \Q2sLob) (solid line) for s = {1, ...,6} and different step sizes h. 
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